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We analyse a specific two dimensional mixed spin Heisenberg model with exchange anisotropy, 
by means of high temperature expansions and Monte Carlo simulations. The goal is to describe 
the magnetic properties of the compound (NBu4)2Mn2 [Cu(opba)]3-6DMSO-H20 which exhibits a 
ferromagnetic transition at Tc — 15K. Extrapolating our analysis on the basis of renormalisation 
group arguments, we find that this transition may result from a very weak anisotropy effect. 

I. INTRODUCTION. 

In the last few years there has been increasing interest in magnetic systems of low dimensionality. For example, 
the rapidly developing field of molecular magnetism deals mainly with quasi one-dimensional and quasi two- 
dimensional compounds. Although the basic theory of their magnetic properties has been known for a long time j|] 
it is now necessary to apply it to the various contexts corresponding to these complex molecular architectures. 

The compound (NBu4)2Mn2[Cu(opba)]3-6DMSO-H20, synthesized by H.O. Stumpf et al [||, exhibits a transition 
at Tc = 15K towards a ferromagnetically ordered state. The structure of this material can be schematically described 
by a superposition of negatively charged layers of hexagonal lattices with the Mn-'^^ ions (spin 5/2) occupying the 
vertices and the Cu''^^ ions (spin 1/2) occupying the middle of the links (see figure 1). The tetrabutylamonium cations, 
NBuJ, are located between the layers. Other compounds of the same kind have been synthesized, differing by the 
nature of the cations between the layers Q. When these cations are small, (Na''", K"*" and tetramethylammonium) a 
long range antiferromagnetic ordering in observed in zero field. An external field of the order of 0.15 kOe is sufficient 
to overcome the very weak interlayer interactions and to lead to a ferromagnetic-like state. The compounds then 
behave as metamagnets. When the cations are larger (tetraethylammonium and beyond), a ferromagnetic ordering 
occurs at a critical temperature Tc- The value of this critical temperature first remains constant and equal to 15K, 
then decreases very smoothly as the cation size increases. In other respects, replacing Mn^^ by a more anisotropic 
spin carrier such as Co^^ results in a significant increase of Tc- These results suggest that both interlayer interactions 
and spin anisotropy are involved in the mechanism of long range ordering. The role of the spin anisotropy in the 
magnetic properties of two-dimensional compounds is much less documented than the three dimensional effects, and 
the goal of this paper is to address this problem. 

In the layer of (NBu4)2Mn2[Cu(opba)]3-6DMSO-Il20, the nearest neighbour Mn-^^ and Cu^^ ions interact through 
an antiferromagnetic coupling. The interlayer interaction in any case is very small as compared to the intralayer one, 
so that, to a good approximation, the spin system can be considered two-dimensional. In a previous paper we have 
shown that such a simple description, in which the 5/2 spins Mn^^ spins are approximated by classical ones, gives a 
good account of the magnetic and thermal properties of the paramagnetic phase of the compound. However, since 
the isotropic 0(3) model is critical only at zero temperature we must include a symmetry breaking mechanism 
in order to explain the phase transition at Tc = 15^. We attribute this symmetry breaking to the presence of spin 
anisotropy. 

Since no single crystal of the Cu-Mn compound has been obtained so far, a direct measurement of the orientation 
of the anisotropy is not possible. However, the existence of a spontaneous magnetization below Tc is the signal of an 
axial anisotropy and of an Ising-like transition. An in-plane anisotropy would have driven the system to be described 
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by an XY symmetric model, which, in two dimensions, exhibits a Kosterhtz-Thouless transition with no ordered phase 

0. 

The aim of this paper, is to investigate the effect of a small axial anisotropy in the simple model described above. 

For weak anisotropy, the critical properties of the model are dominated by the cross-over between the T — critical 
point of the 2D Heisenberg model and the Ising one at Tc- This effect has been widely analysed in the framework 
of the purely classical Heisenberg model |^,^. |l0|jy,|l|,|l|. In particular, the renormalisation group analysis [0 leads 
to the following result : if A is the anisotropy parameter (A = corresponds to the isotropic case) the Ising critical 
temperature decreases to zero as 

T.(A).^ forA^O (1) 

Furthermore, the zero-field susceptibility satisfies a scaling law 

x(A,r)=x(0,T)$(Ae^) (2) 

where the function ^(x) « jcc — for a; « Xc = Ae~. This equation gives additional information on the Ising 

critical region. Let us define the width of this region, S, by 

1-S<^ < l + (5^x(A,T)»x(A = 0,r) 

From eq.(^ we get S w jEAj when A — > 0. Therefore in the weak anisotropy limit we expect the Ising critical region 
to become very narrow and quite close to T = 0. Clearly, this makes experimental investigation difficult. 

On the basis of universality, we transpose these renormalisation group results to our mixed spins system. Since we 
expect the anisotropy to be very small, we need to develop methods specifically designed to handle this cross-over 
effect. In section II we present our techniques of high temperature expansions on the one hand and of Monte-Carlo 
simulations on the other hand. In section III we analyse our results and, from a comparison with the experimental 
data, we determine the value of the anisotropy for the Cu-Mn compound. Conclusions are drawn in the last section. 



II. THE MODEL. 



We denote by sj'^'"^ the spin 



operator associated with the Mn ion at site j, and by S^- " the spin ^ operator 
corresponding to the Cu ion at site i in the middle of a link of the honeycomb lattice. The antiferromagnetic interaction 
is represented by the Heisenberg hamiltonian 



,(Cu) 



i(Mn) 



^^z(Cu)^2(Mn; 
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i=l 



H 
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where J is positive, A(> 0) is the anisotropy parameter, H is the external magnetic field, < i,j > stands for a 
pair of nearest neighbour spins, Ns is the number of sites and is the number of links on the honeycomb lattice 
{Nl — 3/2Ns). The spin | operator can be approximated by a classical spin S.s where s is a unit classical vector 

and S — Id + 1), whereas the spin i operators are expressed in terms of the Pauli matrices, S^^"^ — |cr. Since 
the quantum spin sites are not directly coupled to each other, one can trace out the quantum spin dependence to get 
a completely classical partition function 



Z(r,H)= / I Jl 2 cosh 

\i=i J I <ij> 



/ Ns \ 

exp PgiiisS'R-^Si J 



(4) 



where we have defined 



Wy = --I3JS (s, + + \(sl + 4)2, 



(5) 



and ||X|| stands for the length of the vector X. The indices i and j now label the classical spins located at the vertices 
of the honeycomb lattice. 
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By choosing the orientation of the magnetic field parallel to the anisotropy axis, or orthogonal to it (along the x 
axis), we define Zf^{T, H) = Z{T, He^) with = a; or z, as follows : 



'Ns \ I 1 / Ns 



Z^{T,H)^ J \J[dnA IY12 cosh {K r<^j>) } exp ( PgitisSH ^ j /x = x or z (6) 



<ij> I \ i=l 

with K = yjS and = lis, + s^- + A(sf + - 



<ij> ~ II ' ' 3 ' ^ i ' j'^ JS t^W 

We shall be interested in the standard observables : the specific heat Cy — fc^T^^^ In Z(T, 0), the suceptibility 

, and the total susceptibility, x = 5 Xz + f Xa; i which is measured 

H—O 



along the different directions Xm— In ^ 

experimentally 



III. THE METHOD OF ANALYSIS 



A. The high temperature expansion. 



We have performed the expansion of In in power series of K up to the 19th order and to the second order in H for 
computing the magnetic susceptibility. Then we analysed the series for decreasing values of the anisotropy parameter. 
The complexity of the nearest neighbour interaction (eq.(||)) does not allow the standard techniques HJlJ,!!^ to be 
used. 

The diagrammatic expansion is generated by replacing in eq.(|^) each cosh (^K (j^'^ijy) term by 1 + ^<ij> where the 
function > results from the expansion of the hyperbolic cosine in power series of K to the given maximal order 
and of H to the second order. To each function appearing in the product over the nearest neighbour pairs in 

eq.(|^), is associated one link of a graph G. In this way, a link contains all positive powers of K. 

The powerfuU star graph expansion technique cannot be used here since, due to the presence of the anisotropy 
term, the partition function of articulated graphs does not factorise. Our procedure is based on the the standard 
connected graph expansion for the normalised partition function = Z^(T,H)/ Z^{T,H)\j^q 

\nZ,{T,H) = Y.^{G)u:{G) 

{G} 

where {G} is the set of all connected graphs to a given order, C(G) the embedding constants of the graph G, and 
uj{G) its weight. The weights uj{G) are constructed through the recursive technique \nZ^{G) — X]{g} where {§} 
is the set of all subgraphs of G. 

The main difficulty of the method resides in the computation of InZ^(G), given by : 

lnZ^(G)= / [YldnA (n*n exp(/3.gi/iB5ff5]<J (7) 

where {£} and {1} are respectively the set of links and of vertices belonging to the graph G. We proceed as follows : 

(a) By using the spherical harmonic basis and the reduction formula, each function ^<y> can be expressed 
as : 

iimi/2in2 

where A'^™J is a matrix built recursively of which elements are power series in K and H . 

(b) The exponential term in eq.(j^) is expanded to second order in iJ, and the spin dependence expressed in 
terms of the spherical harmonics. 

(c) The integral is then computed, contracting the products of spherical harmonics by means of the reduction 
formula, and then integrating over the residual angular variables. 
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As an illustration of our results, we give in table I the coefficients of the development of Xz to the zero*'* and first 
order in A. 

The series for the zero-field susceptibility is then analysed by means of the Pade extrapolation technique. For our 
system, the renormalisation group result of eq.(||) becomes 1^ : 

= $(Aetl^) withif=i/3J5 (8) 

where the function $(a;) fa \x — Xc\ ior x fa Xc = Xev^ " . For small A, the singularity is located at large values 
of far away from the perturbative region. According to these results, we analysed the function ^ ^|o'^j ^ in 

terms of the variable w = 1 — e~"^, small at high temperature but bounded for large K . The variational parameter 
a is optimised by stabilising the Pade table. For very low values of A the determination of the critical temperature 
resulting from the location of the poles becomes inaccurate giving large error bars. For A < 10^^ an extension of the 
series to higher orders is needed to reliably get precise results from this method. However, in this very weak regime 
of anisotropy, once the critical temperature is known by another method - Monte-Carlo simulation for instance - the 
series becomes a useful analytic representation of the observables. It will be used to determine the physical parameters 
from a comparison with the experimental data. 

The results of this analysis will be presented and discussed in the next section together with the Monte-Carlo 
results. 



B. The Monte-Carlo simulation. 



We performed a Monte Carlo simulation in order to verify the results of the high temperature expansion, and to 
investigate the very small anisotropy regime where the perturbative technique fails. The simulation is based on the 
effective classical model 

-/3Heff = 51 l^(2cosh||W,j||) (9) 

<i'J> 

where is defined in eq.(|^). The various observables can be expressed as ensemble averages with respect to the 
Boltzmann weight '^"'^ / Z {T , 0). Their expression is a simple generalisation of the definition of eqs.(6,8) of ref. ||] 
in which Wij is replaced by the the new definition of eq.(^. 

Our goal is to explore the weak anisotropy regime where the cross over effect between the 2D Heisenberg and Ising 
fixed points is important. As the anisotopy gets smaller, most of the low temperature region is dominated by the 2D 
Heisenberg regime in which the correlation length remains large due to the essential singularity at T = 0. Therefore, 
the usual limitations of the Monte-Carlo procedure - critical slowing down and the finite size effects - constrain rather 
severely the simulation of these systems. 

In order to overcome the first problem we used a global algorithm. We have adapted the Wolf algorithm to the 
case of an anisotropic interaction. The 'Ising' orientation of the spins, used to contruct the Wolf cluster and which 
is randomly chosen in the standard algorithm, is imposed here by the anisotropy. In our procedure, a Monte-Carlo 
update proceeds in three steps : 

1. construct a first cluster with respect to the z (anisotropy) axis and fiip the corresponding spin components ; 

2. construct another cluster relative to a randomly chosen direction in the x — y plane and flip the corresponding 
spin components ; 

After these two steps, a given spin remains on the cone defined by its initial orientation and the z-axis. Therefore we 
proceed to another step : 

3. change the actual orientation of each spin of the lattice according to a standard Metropolis algorithm. This 
operation is repeated twice. 

We checked that this procedure allows us to recover the results of the zero- anisotopy case ||], and of the strong 
anisotopy (classical Ising) limit with a very small critical slowing down effect. 

The second problem concerns the finite size effects. For small A, since we investigate the low temperature region 
where the correlation length remains large even outside the Ising critical region, we need large lattices. We have 
limited our analysis to A = 0.001, and to a maximum lattice size of L = 256. 
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IV. RESULTS 



A. The critical temperature. 



To validate our methods we compare the resuhs obtained by both procedures. We present in figure 2 the zero-field 
susceptibility, orthogonal to the anisotropy axis (xx) and parallel to it (xz) for the high temperature series and the 
Monte-Carlo data as a function of K and for a moderatedly small value of the anisotropy parameter A = 0.1. The 
agreement between these two results is very good. These curves clearly show the divergent behaviour of the axial 
susceptibility, which will be used to determine the critical temperature in the Monte-Carlo simulation. Figure 3 
displays the specific heat as a function of K for the same anisotropy. Besides the good agreement between the two 
methods, we observe the clear critical signal at Kc — 2.21 which emerges from the comparison with the result for the 
isotropic model. 

The critical temperature is obtained from the Monte-Carlo data by localising the peak in the specific heat and 
the inflexion point of the susceptibility as a function of the temperature. In order to estimate the precision of such 
a determination we performed a complete finite size scaling analysis of the data at A = 0.1 to obtain the critical 
temperature and the susceptibility exponent. The results are presented in figure 4 where wc plotted Tc(L), measured 
from the two signals, for several lattice sizes (fig. (4a)), and in fig. (4b) Inx as a function of In \Kc — K\. Fig. (4a) shows 
a small variation of the finite system critical signal, which allows us to estimate the bulk critical temperature Tc from 
lattices of size not exceeding L = 128. With this estimate, we determine the exponent 7 of the susceptibilty from 
fig. (4b). The result 7 — 1.74(2) is in very good agreement with the expected exact Ising value 7 = 1.75. This is a 
self-consistent indication of the reliability of the critical temperature measurement. Furthermore, we check that the 
Monte-Carlo result falls inside the error bar obtained from the high temperature series analysis. 

With this method, we determined the critical coupling Kc — \ jf^j, as a function of the anisotropy parameter 
on lattices of size L < 256 down to A = 0.001. The results of both analyses - Monte-Carlo and high temperature 
expansion - are displayed in figure 5. Fig. (5a) shows the general trend of the variations of cx as a function of A, 
with a decrease to zero for A ^ 0, and a linear variation at large A. Actually, for A ^ 00, the model coincides with a 
classical Ising model with coupling Rising = so that, in this limit, K~^{X) ~ From the high temperature 

expansion, we get the estimation /^c ising — 1.46(1). The variations for small A are presented in fig. (5b) where we 
plotted Kc vs In A. It appears that the behaviour predicted by the renormalisation group (eq.(|^)) is obtained for very 
small anisotropy. Actually, from fig. (5b) we obtain : 

i^^ = ^ |ln A| + 2.41 A < 0.001. (10a) 
27r 

Therefore, by extrapolating this behaviour down to A — > 0, we are able to predict the critical temperature for very 
small anisotropy. 



B. Comparison with the experimental results. 



In order to compare with the experimental results, we proceed in two steps. 



• We get an estimation of the anisotropy parameter from the value of J previously determined in the analysis 
of the paramagnetic phase |^ and from the experimental critical temperature Tc — 15 K. We obtain Kc — 4.6 
which corresponds to A ~ .0004, according to the fig. (5b). 

• We perform an adjustment of J, 171 and 172 at fixed A by fiting the experimental data with a selected Fade 
approximant. Small variations around the fixed A value do not significantly change the fit. 

The results are presented in figure 6, for the whole range of temperature in fig. (6a) and for the critical region only 
in fig. (6b). An excellent agreement with the experimental data is obtained with the following set of parameters : 

J = 45.5 K , 51 = 2.0 , 32 = 2.11 , A ~ .0005 



The exponential variation of A with respect to Tc (eq.(lOa)) induces rather large error bar on the determination 



of A : actually, an error of one kelvin on Tc induces a variation of A by a factor of 3. However, the magnetic 
susceptibility is rather insensitive to these variations which only affect the Ising critical region, in such a narrow range 
of temperature that it is not visible experimentally. Moreover, the high temperature (T > 50 K) behaviour is left 
completely unchanged by such a small perturbation, thus preserving the good agreement with the isotropic model 
observed in ref.( Q). 
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V. CONCLUSIONS 



We have determined the spin anisotropy which is present in the Cu-Mn magnetic compound |^ and which is 
responsible for a ferromagnetic transition at 15K. Two methods are used in a complementary way, high temperature 
expansion and Monte-Carlo simulation, in order to extract reliable results out of the strong cross-over regime where 
the effect lies. Assuming universality and extending the renormalisation group results to our situation, we obtain a 
very small value for the anisotropy parameter. 

The strong behaviour Tc{X) w is responsible for the fact that a very weak perturbation (A ~ 10~^) produces 
a sizeable effect (Tc ~ 15 K). For these quasi- two dimensional molecular compounds, involving high spin magnetic 
ions, it is unconceivable that at such a low level of magnitude the anisotropy is absent. Therefore we should always 
expect a ferromagnetic transition at an appreciable critical temperature in these systems. 

In this work, we only considered the exchange anisotropy as a source of the 0(3) symmetry breaking. Alternatively, 
the on site anisotropy could be present, but in the limit of weak anisotropy under consideration here, on the ground of 
universality we expect the results to be unchanged. Concerning the interlayer interactions which have been neglected 
in our approach, it is known |]l2| that if the ratio of the interlayer to the intralayer coupling, a = is small, the 

transition temperature induced by these three dimensional effects behaves like Tc « Therefore, even a very 

small interplane coupling may significantly contribute to the measured critical temperature, in competition with the 
anisotropy effect. We have seen that increasing the interlayer distance beyond a certain limit results in a smooth 
decrease of Tc- Taking into account the interlayer coupling would require to introduce an additional cross-over effect 
between the three dimensional model with spatial anisotropy and the two dimensional one with exchange anisotropy 
[ p2[ . As a consequence, the contribution of the spin anisotropy to the critical temperature would be overestimated 
and our result for the anisotropy parameter would turn out to be an upper bound. 
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Table I : A subset of the cocfBcicnts of the high temperature series for the axial susceptibility Xz in powers of 
K = ^PJS up to the 17th order. The coefficient of JsT" is a polynomial of degree n in the anisotropy parameter A. 
For each coefficient we only give the constant and linear term in A (columns with heading A° and A^ respectively) 

The top half of the table contains the even power coefficients ; for each power of A the fraction of the first sub- 
column must be multipied by S^gi'^ and added to the second multiplied by 52^- For instance, the coefficient of K'^ 

is : IS^9i' + W + ^ + W) + O(A^) 

The bottom half of the table contains the odd power coefficients ; the fraction in each column must be multiplied 

by Sgxg2. 
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Figure Captions 



Figure 1 : An elementary cell of the hexagonal lattice, with the Cu ions (closed circles) at the midle of the links and 

the Mn ions (open circles) at the vertices. 

Figure 2 : The susceptibility multiplied by the temperature, as a function of the reduced coupling K = ^fiJS for 
an anisotropy A = 0.1. The data points come from the Monte-Carlo simulation ; the solid line corresponds to 
the Pade approximant of the high temperature expansion (H.T.E.). 

(a) the in-plane susceptibility ; 

(b) the axial susceptibility ; the closed circles correspond to a lattice size L = 128 and the triangles to L = 256. 
Figure 3 : The specific heat as a function of K. 

• for an anisotropy parameter A = 0.1 : the Monte-Carlo data are represented by the closed circles for a lattice 
size L = 128 and by open triangles for L = 256. The solid line corresponds to the high temperature expansion 

(Pade approximant). 

• for the isotropic model (A = 0) : the open circles correspond to the Monte-Carlo simultation [L = 64) and the 
dashed line to the series result.. 

Figure 4 Analysis of the Ising critical regime : 

(a) The finite size critical temperature as a function of the lattice size L : 

• from the maximum of the specific heat (triangles) 

• from the inflexion point of the axial susceptibility (circles) 

(b) the axial susceptibility as as a function of \Kc — K\ in log- log scale. The slope of the fited straight line is 

1.74(2) in good agreement with the expected exact Ising value 7 = 1.75. 

Figure 5 Variation of the critical temperature with A. The series results are represented by open circles and the 

Monte-Carlo simulation by closed circles. 

(a) The large anisotropy regime : we plot the inverse of the reduced critical coupling = ^^jj*^ as a 

function of the anisotropy parameter. The slope of the linear behaviour, expected at large A corresponds 
to incising ^1-46(1). 

(b) The weak anisotropy regime : the reduced critical coupling Kc as a function of A in a logarithmic scale. 

The solid line corresponds to the renormalisation group behaviour Kc = —^\nX + const. The value of 
the constant is 2.41(2). 

Figure 6 Fit of the experimental data (closed circles for sample 1, open circles for sample 2) for the total magnetic 
susceptibility. The solid line corresponds to a Pade approximant of the high temperature series with the 
parameter : J = 45.5 K, gi = 2.0, 32 = 2.14, A = .0005; 

(a) for the whole temperature range 

(b) for the critical region. The difference between the two samples gives an estimation of the experimental errors 
in the critical region. 
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Figure 2 (b) 
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